# Clean up
rm(list = ls())

# Load necessary packages
library(rdd)
library(foreign)
library(tidyverse) 
library(readstata13)
library(gridExtra)
library(psych)
library(readxl)
library(rdrobust)
library(viridis)

# Set working directory: commented; please set your own WD
# setwd("~/Desktop/CPS replication/")

#################################################
######## STUDY 1: COMPARATIVE EVIDENCE ##########  
#################################################

# Import the datasets

RDD <- read.dta13("study1.dta")
df_main_plot <- read_excel("df_main_plot.xlsx", col_names = F)
colnames(df_main_plot) <- c("coefficient", "se", "Model", "dv", "Sample")
df_main_plot$Model <- as.character(df_main_plot$Model)

# Transforming . (missing values in Stata) into NA (missing values in R)

RDD [RDD == "."] <- NA

# Generate a variable for parties above the threshold
RDD$above <- RDD$performance >= 0 

## Analyses ##
# Subset data and remove missing values
rad_right <- RDD[ which(RDD$rad_right == '1'), ]
rad_right <- rad_right[!is.na(rad_right$rad_right),]
rad_right <- rad_right[!is.na(rad_right$performance),]
rad_right <- rad_right[!is.na(rad_right$parl),]

###########
## PLOTS ##
###########

# Descriptive analyses: Showing that radical-right parties are the only ones under-reported

descriptive_data <- RDD %>%
  pivot_longer(c(official_vote, vote_cses), names_to = "source", values_to = "share")

descriptive_matrix <- describeBy(descriptive_data$share, list(descriptive_data$source, descriptive_data$ideology), 
                                 mat = TRUE, digits = 2)

names(descriptive_matrix) [names(descriptive_matrix) == 'group1'] = 'Source'
names(descriptive_matrix) [names(descriptive_matrix) == 'group2'] = 'Ideology'

levels(descriptive_matrix$Source)[levels(descriptive_matrix$Source) == 'vote_cses'] = 'CSES'
levels(descriptive_matrix$Source)[levels(descriptive_matrix$Source) == 'official_vote'] = 'Official'

descriptive_matrix$se = descriptive_matrix$sd/ (sqrt (descriptive_matrix$n))

limits = aes (ymax = mean + (1.96 * se), ymin = mean - (1.96 * se))

descriptive_plot <- descriptive_matrix %>%
  filter(Ideology %in% c("Radical left", "Center left", "Center right", "Radical right")) %>%
  # Generate a variable to order the facets by
  mutate(order = ifelse(Ideology == "Radical left", 1, 
                        ifelse(Ideology == "Center left", 2, 
                               ifelse(Ideology == "Center left", 3,
                                      4)))) %>%
  # Order the facets
  mutate(Ideology = fct_reorder(Ideology, order)) %>%
  ggplot (aes(x = Source, y = mean, fill = Source)) +
  geom_bar (stat = 'identity', position = position_dodge(width = 0.9), alpha = 0.6) +
  geom_errorbar (limits, position = position_dodge(width = 0.9), width = 0) +
  ylab ('Vote share') +
  scale_fill_manual(values = c("blue", "red"), labels = c("Official data", "CSES data")) +
  theme_bw() +
  facet_wrap (~Ideology, nrow = 1, ncol = 4) +
  # Remove labels from x axis
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())
descriptive_plot

ggsave("descriptive_plot.png", plot = descriptive_plot, width = 20, height = 15, units = "cm")

# First stage:
# The discontinuity in the probability of treatment (parliamentary representation), conditional on distance to electoral threshold

first_stage_plot = ggplot(rad_right, aes(x = performance, y = parl, shape = above)) + 
  geom_point(aes(color = as.factor(parl)), alpha = 0.6)  +
  scale_color_manual(values = c("red","blue")) +
  ylim(0, 1) + 
  xlim(-5, 5) + 
  stat_smooth(method = "loess",span = 1, colour = "black", size = 0.75, se = F) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") +
  xlab("Distance to the threshold (percentage points)") +
  ylab("Parliamentary representation") + 
  theme(legend.position = "none", panel.background = element_rect(fill = "white", colour = "grey50"))
first_stage_plot

ggsave("firststage.png", plot = first_stage_plot, width = 20, height = 15, units = "cm")

# Main plot: discontinuity in the proportion of the official vote reported in CSES, conditional on distance to threshold
rd_plot <- ggplot(rad_right, aes(x = performance, y = reported_vote)) +
  scale_x_continuous(limits = c(-5, 5)) +
  coord_cartesian(ylim = c(-0.1, 2)) +
  stat_summary_bin( aes(colour = as.factor(above)), fun.y = 'mean', bins = 20,
                    color = 'black', size = 0.75, geom = 'point', alpha = 0.6) +
  stat_smooth(method = "loess", aes(colour = as.factor(above)), span = 1.5, size = 0.5) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") +
  xlab("Difference between vote share and threshold (percentage points)") +
  ylab("Proportion of the official vote for RRPs reported in post-electoral surveys") + 
  scale_colour_manual(name = " ", values = c("blue", "blue")) +
  theme(legend.position = "none", 
        panel.background = element_rect(fill = "white", colour = "grey50")) 
rd_plot

# Plot the coefficients
main_plot_coefs <- df_main_plot %>%
  mutate(ci_upper_95 = coefficient + 1.96 * se) %>%
  mutate(ci_lower_95 = coefficient - 1.96 * se) %>%  
  mutate(ci_upper_90 = coefficient + 1.645 * se) %>%
  mutate(ci_lower_90 = coefficient - 1.645 * se) %>%
  mutate(order = ifelse(Sample == "Radical right", 0, 1)) %>%
  mutate(Sample = fct_reorder(Sample, order)) %>%
  ggplot(aes(Model, coefficient, color = dv, shape = dv)) + 
  geom_hline(aes(yintercept = 0), linetype = 2, color = "gray20") + 
  coord_flip(ylim = c(-1.5, 1.5)) + 
  geom_linerange(aes(ymin = ci_lower_95, ymax = ci_upper_95), position = position_dodge(width = 0.75), size = 0.5) + 
  geom_linerange(aes(ymin = ci_lower_90, ymax = ci_upper_90), position = position_dodge(width = 0.75), size = 1) + 
  geom_point(size = 1.75, position = position_dodge(width = 0.75)) + 
  theme_minimal( ) +
  theme(legend.title = element_blank(),
        legend.position = c(0.8, 0.05),
        panel.background = element_rect(fill = NA, color = "black")) +
  ylab("LATE") + 
  xlab("Model") +
  facet_grid(Sample ~ ., scale = 'free') +
  scale_color_manual(values = c("red", "blue")) +
  guides(colour = guide_legend(reverse = TRUE), shape = guide_legend(reverse = TRUE))
main_plot_coefs

main_plot_final <- grid.arrange(
  rd_plot,
  main_plot_coefs,
  nrow = 1,
  ncol = 2)
main_plot_final

ggsave("main_plot.png", plot = main_plot_final, width = 30, height = 15, units = "cm")

# RDD plot for center right parties
cr_plot <- RDD %>%
  filter(ches_genlr >= 5.5 & ches_genlr <= 8 & rad_right == 0) %>%
  ggplot(aes(x = performance, y = reported_vote)) +
  scale_x_continuous(limits = c(-5, 5)) +
  coord_cartesian(ylim = c(-0, 2.5)) +
  stat_summary_bin(aes(colour = as.factor(above)), fun.y = 'mean', bins = 20,
                   color = 'black', size = 0.75, geom = 'point', alpha = 0.6) +
  stat_smooth(method = "loess", aes(colour = as.factor(above)), span = 1.5, size = 0.5) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") +
  xlab("Difference between vote share and threshold (percentage points)") +
  ylab("Proportion of the official vote for RRPs reported in post-electoral surveys") + 
  scale_colour_manual(name = " ", values = c("blue", "blue")) +
  theme(legend.position = "none", 
        panel.background = element_rect(fill = "white", colour = "grey50")) 
cr_plot

ggsave("cr_plot.png", plot = cr_plot, width = 30, height = 20, units = "cm")

## CHECKING WHO MOVES ##
# Main RD with reported vote for center right parties as outcome
movers_cr <- rad_right %>%
  filter(max_dummy == 1) %>%
  ggplot(aes(x = performance, y = reported_vote_total_cr)) +
  scale_x_continuous(limits = c(-5, 5)) +
  coord_cartesian(ylim = c(-0, 2.5)) +
  stat_summary_bin(aes(colour = as.factor(above)), fun.y = 'mean', bins = 20,
                   color = 'black', size = 0.75, geom = 'point', alpha = 0.6) +
  stat_smooth(method = "loess", aes(colour = as.factor(above)), span = 1.5, size = 0.5) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") +
  xlab("Difference between vote share for RRP and threshold (percentage points)") +
  ylab("Proportion of the official vote reported in post-electoral surveys") + 
  scale_colour_manual(name = " ", values = c("blue", "blue")) +
  theme(legend.position = "none", 
        panel.background = element_rect(fill = "white", colour = "grey50")) +
  ggtitle ("Center-right parties")
movers_cr

ggsave("cr_plot.png", plot = cr_plot, width = 30, height = 20, units = "cm")

# Main RD with reported vote for center left parties as outcome
movers_cl <- rad_right %>%
  filter(max_dummy == 1) %>%
  ggplot(aes(x = performance, y = reported_vote_total_cl)) +
  scale_x_continuous(limits = c(-5, 5)) +
  coord_cartesian(ylim = c(-0, 2.5)) +
  stat_summary_bin(aes(colour = as.factor(above)), fun.y = 'mean', bins = 20,
                   color = 'black', size = 0.75, geom = 'point', alpha = 0.6) +
  stat_smooth(method = "loess", aes(colour = as.factor(above)), span = 1.5, size = 0.5) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") +
  xlab("Difference between vote share for RRP and threshold (percentage points)") +
  ylab("Proportion of the official vote reported in post-electoral surveys") + 
  scale_colour_manual(name = " ", values = c("red", "red")) +
  theme(legend.position = "none", 
        panel.background = element_rect(fill = "white", colour = "grey50")) +
  ggtitle ("Center-left parties")
movers_cl

# Main RD with reported vote for all left parties as outcome
movers_l <- rad_right %>%
  filter(max_dummy == 1) %>%
  ggplot(aes(x = performance, y = reported_vote_total_l)) +
  scale_x_continuous(limits = c(-5, 5)) +
  coord_cartesian(ylim = c(-0, 2.5)) +
  stat_summary_bin(aes(colour = as.factor(above)), fun.y = 'mean', bins = 20,
                   color = 'black', size = 0.75, geom = 'point', alpha = 0.6) +
  stat_smooth(method = "loess", aes(colour = as.factor(above)), span = 1.5, size = 0.5) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") +
  xlab("Difference between vote share for RRP and threshold (percentage points)") +
  ylab("Proportion of the official vote reported in post-electoral surveys") + 
  scale_colour_manual(name = " ", values = c("orange", "orange")) +
  theme(legend.position = "none", 
        panel.background = element_rect(fill = "white", colour = "grey50")) +
  ggtitle ("All left parties")
movers_l

# Main RD with reported turnout
movers_turnout <- rad_right %>% 
  filter(max_dummy == 1) %>%
  ggplot(aes(x = performance, y = reported_vote_turnout)) +
  scale_x_continuous(limits = c(-5, 5)) +
  coord_cartesian(ylim = c(-0, 2.5)) +
  stat_summary_bin(aes(colour = as.factor(above)), fun.y = 'mean', bins = 20,
                   color = 'black', size = 0.75, geom = 'point', alpha = 0.6) +
  stat_smooth(method = "loess", aes(colour = as.factor(above)), span = 1.5, size = 0.5) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") +
  xlab("Difference between vote share for RRP and threshold (percentage points)") +
  ylab("Proportion of the official turnout reported in post-electoral surveys") + 
  scale_colour_manual(name = " ", values = c("black", "black")) +
  theme(legend.position = "none", 
        panel.background = element_rect(fill = "white", colour = "grey50")) +
  ggtitle ("Turnout")
movers_turnout

# Main RD with reported turnout
movers_dkna <- rad_right %>% 
  filter(max_dummy == 1) %>%
  ggplot(aes(x = performance, y = dkna)) +
  scale_x_continuous(limits = c(-5, 5)) +
  #  coord_cartesian(ylim = c(-0, 2.5)) +
  stat_summary_bin(aes(colour = as.factor(above)), fun.y = 'mean', bins = 20,
                   color = 'black', size = 0.75, geom = 'point', alpha = 0.6) +
  stat_smooth(method = "loess", aes(colour = as.factor(above)), span = 1.5, size = 0.5) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") +
  xlab("Difference between vote share for RRP and threshold (percentage points)") +
  ylab("Proportion of the respondents not answering who they voted for") + 
  scale_colour_manual(name = " ", values = c("darkgreen", "darkgreen")) +
  theme(legend.position = "none", 
        panel.background = element_rect(fill = "white", colour = "grey50")) +
  ggtitle ("Don't know / Don't Answer")
movers_dkna

# Putting the three plots together
movers_all <- grid.arrange(
  movers_cr,
  movers_cl,
  movers_l,
  movers_turnout,
  movers_dkna,
  nrow = 2,
  ncol = 3)
movers_all

ggsave("rdmovers.png", plot = movers_all, width = 45, height = 30, units = "cm")

# Analyses with country FE
rrbw <- IKbandwidth(rad_right$performance, rad_right$reported_vote, verbose=TRUE)

rd_rr <- RDestimate(reported_vote ~ performance + parl | cntr1 + cntr2 + cntr3 + cntr4 + 
                      cntr5 + cntr6 + cntr7 + cntr8 + cntr9 + cntr10 + cntr11 + cntr12 + 
                      cntr14 + cntr15 + cntr16 + cntr17 + cntr18 + cntr20, 
                    cluster = rad_right$country_code, 
                    data = rad_right,
                    bw = rrbw) 
summary(rd_rr)

# Looking for sorting around the threshold
sorting <- rad_right[rad_right$performance <= 5 & rad_right$performance >= -5 ,]
DCdensity(sorting$performance, plot = TRUE, bin = 0.5)
title(xlab = "Centered Vote Share")
abline(v = 0, lty = 1, col = "grey50")

#################################################
## Analyses depending on size of the threshold ##
#################################################

high_thresholds <- rad_right %>%
  filter(high_threshold == 1) %>%
  ggplot(aes(x = performance, y = reported_vote)) +
  scale_x_continuous(limits = c(-5, 5)) +
  coord_cartesian(ylim = c(-0.1, 2)) +
  stat_summary_bin( aes(colour = as.factor(above)), fun.y = 'mean', bins = 20,
                    color = 'black', size = 0.75, geom = 'point', alpha = 0.6) +
  stat_smooth(method = "loess", aes(colour = as.factor(above)), span = 1.5, size = 0.5) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") +
  xlab("Difference between vote share for an RRP and threshold (percentage points)") +
  ylab("Proportion of the official vote for other RRP reported in post-electoral surveys") + 
  scale_colour_manual(name = " ", values = c("blue", "blue")) +
  ggtitle("High thresholds") +
  theme(legend.position = "none", 
        panel.background = element_rect(fill = "white", colour = "grey50")) 
high_thresholds

low_thresholds <- rad_right %>%
  filter(high_threshold == 0) %>%
  ggplot(aes(x = performance, y = reported_vote)) +
  scale_x_continuous(limits = c(-5, 5)) +
  coord_cartesian(ylim = c(-0.1, 2)) +
  stat_summary_bin( aes(colour = as.factor(above)), fun.y = 'mean', bins = 20,
                    color = 'black', size = 0.75, geom = 'point', alpha = 0.6) +
  stat_smooth(method = "loess", aes(colour = as.factor(above)), span = 1.5, size = 0.5) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") +
  xlab("Difference between vote share for an RRP and threshold (percentage points)") +
  ylab("Proportion of the official vote for other RRP reported in post-electoral surveys") + 
  ggtitle("Low thresholds") +
  scale_colour_manual(name = " ", values = c("red", "red")) +
  theme(legend.position = "none", 
        panel.background = element_rect(fill = "white", colour = "grey50")) 
low_thresholds

# Putting together the two plots
thresholds <- grid.arrange(
  low_thresholds,
  high_thresholds,
  nrow = 1,
  ncol = 2)
thresholds

ggsave("thresholds.png", plot = thresholds, width = 30, height = 15, units = "cm")

##################################################################################
## Descriptive analyses depending on whether the party was in parliament before ##
##################################################################################

previously_in_parl <- rad_right %>%
  filter(inparl_before == 1) %>%
  ggplot(aes(x = performance, y = reported_vote)) +
  scale_x_continuous(limits = c(-5, 5)) +
  coord_cartesian(ylim = c(-0.1, 2)) +
  stat_summary_bin( aes(colour = as.factor(above)), fun.y = 'mean', bins = 20,
                    color = 'black', size = 0.75, geom = 'point', alpha = 0.6) +
  stat_smooth(method = "loess", aes(colour = as.factor(above)), span = 1.5, size = 0.5) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") +
  xlab("Difference between vote share for an RRP and threshold (percentage points)") +
  ylab("Proportion of the official vote for other RRP reported in post-electoral surveys") + 
  scale_colour_manual(name = " ", values = c("red", "red")) +
  ggtitle("Radical right parties in parliament in the previous election") +
  theme(legend.position = "none", 
        panel.background = element_rect(fill = "white", colour = "grey50")) 
previously_in_parl

previously_not_in_parl <- rad_right %>%
  filter(inparl_before == 0) %>%
  ggplot(aes(x = performance, y = reported_vote)) +
  scale_x_continuous(limits = c(-5, 5)) +
  coord_cartesian(ylim = c(-0.1, 2)) +
  stat_summary_bin( aes(colour = as.factor(above)), fun.y = 'mean', bins = 20,
                    color = 'black', size = 0.75, geom = 'point', alpha = 0.6) +
  stat_smooth(method = "loess", aes(colour = as.factor(above)), span = 1.5, size = 0.5) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") +
  xlab("Difference between vote share for an RRP and threshold (percentage points)") +
  ylab("Proportion of the official vote for other RRP reported in post-electoral surveys") + 
  ggtitle("Radical right parties not in parliament in the previous election") +
  scale_colour_manual(name = " ", values = c("blue", "blue")) +
  theme(legend.position = "none", 
        panel.background = element_rect(fill = "white", colour = "grey50")) 
previously_not_in_parl

# Putting together the two plots
previoustreat <- grid.arrange(
  previously_in_parl,
  previously_not_in_parl,
  nrow = 1,
  ncol = 2)
previoustreat

ggsave("previoustreat.png", plot = previoustreat, width = 30, height = 15, units = "cm")

#####################################
## Looking at effects through time ##
#####################################

longitudinal_both <- read.dta("predictboth.dta")

longitudinal_both$parl <- as.factor(longitudinal_parl$parl)

longitudinal_both_plot <- longitudinal_both %>%
  ggplot(aes(x = dse, y = predvalue, color = as.factor(parl))) + 
  geom_point(size = 1.4) +
  scale_x_continuous(name = "Number of days between election and survey interview",limits = c(0, 80)) +
  scale_y_continuous(name = "Predicted probability of reporting having voted for a RRP")+ 
  geom_hline(yintercept = 0, 
             linetype = 2, color = "red") + 
  guides(fill = guide_legend(reverse = TRUE), colour = guide_legend(reverse = TRUE)) +
  ggtitle(" ") +
  scale_color_manual(name = " ",labels = c("RRP out of parliament", "RRP in parliament"), values = c("black","red")) + 
  stat_smooth(method = "loess", formula = y ~ x, size = 0.5, span = 1) +
  theme_minimal(base_size = 14)  + 
  theme(legend.position = c(0.15, .9), legend.key.size = unit(1, "line"),
        axis.text = element_text(size = 14),axis.title.x = element_text(size = 14),axis.title.y = element_text(size = 14),
        plot.title = element_text(size = 14, hjust = 0.5))
longitudinal_both_plot

ggsave("longitudinal.png", plot = longitudinal_both_plot, width = 20, height = 15, units = "cm")

## Checking if there is evolution in the proportion of official vote that is reported in post-electoral surveys across time ##
repvote <- read.dta13("repvote_pred.dta")

repvote_plot <- ggplot(data = repvote, aes(x = year, y = repvote_pred)) + 
  geom_point(size = 1.4) +
  theme(legend.position = "none") +
  ylab("Proportion of the vote for RRPs that is reported in post-electoral surveys")  + 
  scale_x_continuous(name="Year",limits = c(2001, 2018)) +
  scale_y_continuous(limits = c(0, 1.5)) +
  theme(legend.position = "none") +    
  stat_smooth(method = "loess", formula = y ~ x, size = 0.5, span = 1, color = "black") +
  geom_hline(yintercept = 0, 
             linetype = 2, color = "red") +
  theme_bw()
repvote_plot  

ggsave("repvote_time.png", plot = repvote_plot,  width = 20, height = 15, units = "cm") 

########################
# Jackknife estimates #
########################

# Import dataset with results
jackknives <- read_excel("jackknives.xlsx", col_names = F)

# Add column names
colnames(jackknives) <- c("iteration number", "coefficient", "se", "model", "order")

# Add a column with the country that has been excluded in each model
country <- rep(c("Excluding Austria", "Excluding Bulgaria", "Excluding Czech Rep", 
                 "Excluding Germany", "Excluding Denmark", "Excluding Estonia", 
                 "Excluding Greece", "Excluding Croatia", "Excluding Hungary", 
                 "Excluding Israel", "Excluding Italy", "Excluding Latvia", 
                 "Excluding Netherlands", "Excluding Norway", "Excluding Poland",
                 "Excluding Romania", "Excluding Serbia", "Excluding Slovakia", 
                 "Excluding Slovenia", "Excluding Sweden", "Excluding Ukraine"),
               times = 4)
jk <- cbind(jackknives, country)

# set up theme for plots
theme_base1 <- 
  theme_minimal(base_size=9)  + 
  theme(legend.position = c(0.4, 0.2),
        legend.direction = "horizontal",
        legend.key.size = unit(0.5, "cm"), 
        axis.text=element_text(size=12),axis.title.x=element_text(size=12),axis.title.y=element_text(size=10),
        plot.title = element_text(size=12, hjust= 0.5))

# Make the actual plot
jk_plot <- jk %>%
  mutate(model = fct_reorder(as.factor(model), order)) %>%
  mutate(ci_upper95 = coefficient + 1.96 * se) %>%
  mutate(ci_lower95 = coefficient - 1.96 * se) %>%
  mutate(ci_upper90 = coefficient + 1.645 * se) %>%
  mutate(ci_lower90 = coefficient - 1.645 * se) %>%
  ggplot(aes(x = country, y = coefficient)) +
  geom_linerange(aes(ymin = ci_lower95, ymax = ci_upper95), position = position_dodge(width = 0.75), size = 0.5) + 
  geom_linerange(aes(ymin = ci_lower90, ymax = ci_upper90), position = position_dodge(width = 0.75), size = 1) + 
  geom_point(size = 1.25, position = position_dodge(width = 0.75) ) + 
  geom_hline(aes(yintercept = 0), linetype = 2, color = "red") + 
  coord_flip() + 
  theme_bw() +
  xlab(" ") + 
  ylab("LATE") + 
  facet_grid(~model) +
  ggtitle(" ")
jk_plot

ggsave("jk.png", plot = jk_plot, width = 35, height = 18, units = "cm") 

###########################################
# Changing LATE with Different Bandwidths #
###########################################

bw <- seq (2.5, 7.5, 0.1)
est_conventional <- rep (NA, length(bw))
se_conventional <- rep (NA, length(bw))
rdd.list_conventional <- list ()

for(i in 1:length(bw)){
  rdd.list_conventional[[i]] <- rdrobust(y = rad_right$reported_vote, x = rad_right$performance, 
                                         h = bw[i], 
                                         fuzzy = rad_right$parl, 
                                         cluster = rad_right$country_code)
  est_conventional[i] <- rdd.list_conventional[[i]]$coef[1]
  se_conventional[i] <- rdd.list_conventional[[i]]$se[1]
}

df_conventional <- data.frame(
  bw = bw,
  est = est_conventional,
  ci.hi95 = est_conventional + 1.96*se_conventional,
  ci.lo95 = est_conventional - 1.96*se_conventional,
  ci.hi90 = est_conventional + 1.645*se_conventional,
  ci.lo90 = est_conventional - 1.645*se_conventional)

bw_conventional_plot <- ggplot(df_conventional, aes(x = bw, y = est)) +  
  geom_line(aes(x = bw, y = est_conventional), size = 1, color = "blue") + 
  geom_line(aes(x = bw, y = ci.hi90), color = "black") +
  geom_line(aes(x = bw, y = ci.lo90), color = "black") +
  geom_line(aes(x = bw, y = ci.hi95), linetype = "dashed", color = "black") +
  geom_line(aes(x = bw, y = ci.lo95), linetype = "dashed", color = "black") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  #  geom_vline(xintercept = 3.918, linetype = "dashed",color = "green4",size = 1) +
  xlab("Bandwidth") +
  ggtitle("Conventional procedure") +
  ylab("LATE")   + 
  theme_bw() +
  scale_x_continuous(breaks = c(2,3,4,5,6,7,8)) + 
  scale_y_continuous(breaks = c(-0.3,-0.2,-0.1,0,0.1, 0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1, 1.1, 1.2), limits = c(-0.3,1.2)) + 
  theme(plot.title = element_text(face="bold", size = 13))
bw_conventional_plot

est_biascorrect <- rep (NA, length(bw))
se_biascorrect <- rep (NA, length(bw))
rdd.list_biascorrect <- list ()

for(i in 1:length(bw)){
  rdd.list_biascorrect[[i]] <- rdrobust(y = rad_right$reported_vote, x = rad_right$performance, 
                                        h = bw[i], 
                                        fuzzy = rad_right$parl, 
                                        cluster = rad_right$country_code)
  est_biascorrect[i] <- rdd.list_biascorrect[[i]]$coef[2]
  se_biascorrect[i] <- rdd.list_biascorrect[[i]]$se[2]
}

df_biascorrect <- data.frame(
  bw = bw,
  est = est_biascorrect,
  ci.hi95 = est_biascorrect + 1.96*se_biascorrect,
  ci.lo95 = est_biascorrect - 1.96*se_biascorrect,
  ci.hi90 = est_biascorrect + 1.645*se_biascorrect,
  ci.lo90 = est_biascorrect - 1.645*se_biascorrect)

bw_biascorrect_plot <- ggplot(df_biascorrect, aes(x = bw, y = est)) +  
  geom_line(aes(x = bw, y = est), size = 1, color = "blue") + 
  geom_line(aes(x = bw, y = ci.hi90), color = "black") +
  geom_line(aes(x = bw, y = ci.lo90), color = "black") +
  geom_line(aes(x = bw, y = ci.hi95), linetype = "dashed", color = "black") +
  geom_line(aes(x = bw, y = ci.lo95), linetype = "dashed", color = "black") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  #  geom_vline(xintercept = 3.918, linetype = "dashed",color = "green4",size = 1) +
  xlab("Bandwidth") +
  ggtitle("Bias-corrected procedure") +
  ylab("LATE")   + 
  theme_bw() +
  scale_x_continuous(breaks = c(2,3,4,5,6,7,8)) + 
  scale_y_continuous(breaks = c(-0.3,-0.2,-0.1,0,0.1, 0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1, 1.1, 1.2), limits = c(-0.3,1.2)) + 
  theme(plot.title = element_text(face="bold", size = 13))
bw_biascorrect_plot

est_robust <- rep (NA, length(bw))
se_robust <- rep (NA, length(bw))
rdd.list_robust <- list ()

for(i in 1:length(bw)){
  rdd.list_robust[[i]] <- rdrobust(y = rad_right$reported_vote, x = rad_right$performance, 
                                   h = bw[i], 
                                   fuzzy = rad_right$parl, 
                                   cluster = rad_right$country_code)
  est_robust[i] <- rdd.list_robust[[i]]$coef[3]
  se_robust[i] <- rdd.list_robust[[i]]$se[3]
}

df_robust <- data.frame(
  bw = bw,
  est = est_robust,
  ci.hi95 = est_robust + 1.96*se_robust,
  ci.lo95 = est_robust - 1.96*se_robust,
  ci.hi90 = est_robust + 1.645*se_robust,
  ci.lo90 = est_robust - 1.645*se_robust)

bw_robust_plot <- ggplot(df_robust, aes(x = bw, y = est)) +  
  geom_line(aes(x = bw, y = est), size = 1, color = "blue") + 
  geom_line(aes(x = bw, y = ci.hi90), color = "black") +
  geom_line(aes(x = bw, y = ci.lo90), color = "black") +
  geom_line(aes(x = bw, y = ci.hi95), linetype = "dashed", color = "black") +
  geom_line(aes(x = bw, y = ci.lo95), linetype = "dashed", color = "black") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  #  geom_vline(xintercept = 3.918, linetype = "dashed",color = "green4",size = 1) +
  xlab("Bandwidth") +
  ggtitle("Robust procedure") +
  ylab("LATE")   + 
  theme_bw() +
  scale_x_continuous(breaks = c(2,3,4,5,6,7,8)) + 
  scale_y_continuous(breaks = c(-0.3,-0.2,-0.1,0,0.1, 0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1, 1.1, 1.2), limits = c(-0.3,1.2)) + 
  theme(plot.title = element_text(face="bold", size = 13))
bw_robust_plot

# Putting together the three plots
rdbw_all <- grid.arrange(
  bw_conventional_plot,
  bw_biascorrect_plot,
  bw_robust_plot,
  nrow = 1,
  ncol = 3)
rdbw_all

ggsave("rdbw.png", plot = rdbw_all, width = 25, height = 10, units = "cm")

## Checking if there is association between left-right mean placement of remaining parties and reported vote variable ##

lrmean <- ggplot(data = rad_right, aes(x = lrmean_party, y = reported_vote)) + 
  geom_point(size = 1.4) +
  xlab("Mean left-right placement of non-radical right parties") + theme(legend.position = "none") +
  ylab("Proportion of the official vote that is reported in post electoral surveys")  + theme(legend.position = "none") +
  geom_hline(yintercept = 0, 
             linetype = 2, color = "red") +
  geom_smooth(method = "loess", formula = y ~ x, size = 0.5, span = 1, color = "black") +
  theme_base1
lrmean

ggsave("lrmean.png", plot = lrmean,  width = 15, height = 12.5, units = "cm")  

###############################################################
####### STUDY 2: COMPARING DIFFERENT MODES OF INTERVIEW ####### 
###############################################################

# Import data
ppfull <- read.dta13("study2.dta")

# Import dataset
study2 <- read_excel("study2.xlsx", col_names = F)

# Give columns the correct names
colnames(study2) <- c("group", "coef", "se", "Parliament", "rad_right", "Interview")

study2$Interview <- as.factor(study2$Interview)

# Fix level values
levels(study2$Interview)[levels(study2$Interview) == '1'] = 'Public'
levels(study2$Interview)[levels(study2$Interview) == '0'] = 'Private'

levels(study2$Parliament)[levels(study2$Parliament) == '0'] = 'No'
levels(study2$Parliament)[levels(study2$Parliament) == '1'] = 'Yes'

levels(study2$rad_right)[levels(study2$rad_right) == '0'] = 'Remaining Parties'
levels(study2$rad_right)[levels(study2$rad_right) == '1'] = 'Radical Right Parties'

# Make the plot
interview <- study2 %>%
  mutate(order = ifelse(Interview == "Private", 1, 0)) %>%
  mutate(Interview = fct_reorder(Interview, order)) %>%
  mutate(parl = ifelse(Parliament == 1, "Yes", "No")) %>%
  mutate(rrp = ifelse(rad_right == 1, "Radical right parties", "Remaining parties")) %>%
  mutate(ci_upper = coef + 1.96 * se) %>%
  mutate(ci_lower = coef - 1.96 * se) %>%
  ggplot (aes(x = parl, y = coef, fill = Interview)) +
  geom_bar (stat = 'identity',  position = position_dodge(0.9), alpha = 0.6) +
  geom_errorbar(aes(ymin = ci_lower, ymax = ci_upper), width = 0, position = position_dodge(0.9)) +  
  labs (x = 'Parliamentary representation',
        y ='Reported like / dislike for the Party') +
  theme_bw() +
  scale_fill_viridis(discrete = TRUE) +
  facet_wrap(~rrp)
interview

ggsave("interview.png", plot = interview, width = 20, height = 15, units = "cm")

#################################################################
## Looking for heterogenity based on vote in previous election ##
#################################################################

# Import dataset
study2_left <- read_excel("study2_left.xlsx", col_names = F)

# Give columns the correct names
colnames(study2_left) <- c("group", "coef", "se", "parl", "left_voter", "Interview")

study2_left$Interview <- as.factor(study2_left$Interview)

# Fix level values
levels(study2_left$Interview)[levels(study2_left$Interview) == '1'] = 'Public'
levels(study2_left$Interview)[levels(study2_left$Interview) == '0'] = 'Private'

# Make the plot
study2_left_plot <- study2_left %>%
  mutate(order = ifelse(Interview == "Private", 1, 0)) %>%
  mutate(Interview = fct_reorder(Interview, order)) %>%
  mutate(Parliament = ifelse(parl == 1, "Yes", "No"))  %>%
  mutate(left = ifelse(left_voter == 1, "Left voters", "Right voters"))  %>%
  mutate(ci_upper = coef + 1.96 * se) %>%
  mutate(ci_lower = coef - 1.96 * se) %>%
  mutate(order = ifelse(Parliament == "RRP in parliament", 1, 0)) %>%
  mutate(Parliament = fct_reorder(Parliament, order)) %>%
  ggplot (aes(x = Parliament, y = coef, fill = Interview)) +
  geom_bar (stat = 'identity', position = position_dodge(0.9), alpha = 0.6) +
  geom_errorbar(aes(ymin = ci_lower, ymax = ci_upper), width = 0, position = position_dodge(0.9)) +  
  labs (x = 'Parliamentary representation',
        y ='Reported like / dislike for the RRP') +
  theme_bw() +
  scale_fill_viridis(discrete = TRUE) +
  facet_wrap(~left)
study2_left_plot

ggsave("movers.png", plot = study2_left_plot, width = 20, height = 15, units = "cm")

#######
# Analyses of heterogeneity based on college degree
#######

# Import dataset
study2_college <- read_excel("study2_college.xlsx", col_names = F)

# Give columns the correct names
colnames(study2_college) <- c("group", "coef", "se", "parl", "college", "Interview")

study2_college$Interview <- as.factor(study2_college$Interview)

# Fix level values
levels(study2_college$Interview)[levels(study2_college$Interview) == '1'] = 'Public'
levels(study2_college$Interview)[levels(study2_college$Interview) == '0'] = 'Private'

levels(study2_college$parl)[levels(study2_college$parl) == '0'] = 'No'
levels(study2_college$parl)[levels(study2_college$parl) == '1'] = 'Yes'

levels(study2_college$college)[levels(study2_college$college) == '0'] = 'No'
levels(study2_college$college)[levels(study2_college$college) == '1'] = 'Yes'

# Make the plot
study2_college_plot <- study2_college %>%
  mutate(order = ifelse(Interview == "Private", 1, 0)) %>%
  mutate(Interview = fct_reorder(Interview, order)) %>%
  mutate(College = ifelse(college == 1, "Yes", "No")) %>%
  mutate(Parliament = ifelse(parl == 1, "RRP in parliament", "RRP out of parliament"))  %>%
  mutate(ci_upper = coef + 1.96 * se) %>%
  mutate(ci_lower = coef - 1.96 * se) %>%
  mutate(order = ifelse(Parliament == "RRP in parliament", 1, 0)) %>%
  mutate(Parliament = fct_reorder(Parliament, order)) %>%
  ggplot (aes(x = College, y = coef, fill = Interview)) +
  geom_bar (stat = 'identity', position = position_dodge(0.9), alpha = 0.6) +
  geom_errorbar(aes(ymin = ci_lower, ymax = ci_upper), width = 0, position = position_dodge(0.9)) +  
  labs (x = 'College Degree',
        y ='Reported like / dislike for the RRP') +
  theme_bw() +
  scale_fill_viridis(discrete = TRUE) +
  facet_wrap(~Parliament)
study2_college_plot

ggsave("het_col.png", plot = study2_college_plot, width = 20, height = 15, units = "cm")

###########################################
####### STUDY 3: CASE STUDY OF UKIP ####### 
###########################################

# Parallel trends graph
# Importing dataset
parallel <- read.dta13("partrends_genelec.dta", encoding = 'UTF-8')

# Making sure all numeric variables are such
parallel$meanrepvote <- as.numeric(parallel$meanrepvote)
parallel$evertreat <- as.factor(parallel$evertreat)
parallel$time_period <- as.numeric(parallel$time_period)

cbPalette <- c("red","black")

partrends = ggplot(parallel, aes(x = time_period, y = meanrepvote, colour = evertreat)) +  
  geom_line(aes(linetype = evertreat, color = evertreat), position = position_dodge(0.25)) + 
  geom_point(position = position_dodge(0.25)) +
  scale_y_continuous(limits = c(0, 1.3), breaks = c(0, 0.25, 0.5, 0.75, 1, 1.25))  + 
  ylab("Proportion of the official vote reported in the BES") +
  scale_color_manual(name = "Type of Party",values = cbPalette,labels = c( ), guide = "none") +
  scale_linetype_manual(name = " ", values = c(1,2), guide = "none")+
  scale_x_continuous(name = "Election", breaks = c(1, 2, 3), 
                     labels = c("2005", "2010", "2015")) +
  annotate("text", x = 2, y = 1.16, label = "Control: Conservative, Labour and LibDem", size = 3) +
  annotate("text", x = 2, y = 0.58, label = "Treated: UKIP", size = 3) +
  labs(fill = " ") +
  theme_classic() + 
  geom_errorbar(aes(ymin = ci_lower, ymax = ci_upper), width = 0.1, position = position_dodge(0.25)) +  
  geom_vline(xintercept = 2.85, linetype = "dashed",color = "grey",size = .5)
partrends

# Plotting the findings of the DiD model

did <- read.dta13("ukip_graph.dta", encoding = 'UTF-8')

# Adding line breaks
levels(did$Model) <- gsub("just pre", "just\npre", levels(did$Model))
levels(did$Model) <- gsub("With demographic", "With\ndemographic", levels(did$Model))
levels(did$Model) <- gsub("demographic weights", "demographic\nweights", levels(did$Model))

did_plot <- did %>%
  filter (SE == "Clustered") %>%
  mutate(Model = fct_reorder(Model, desc(order))) %>%
  ggplot (aes(x = Model, y = coef)) +
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +
  geom_linerange(aes(ymin = ci_lower_95, ymax = ci_upper_95), position = position_dodge(width = 0.75), size = 0.5) + 
  geom_linerange(aes(ymin = ci_lower_90, ymax = ci_upper_90), position = position_dodge(width = 0.75), size = 1) + 
  geom_point(size = 1, position = position_dodge(width = 0.75) ) + 
  scale_y_continuous(name = "ATT", limits = c(-0.6, 0.6), breaks = c(-0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6)) +
  coord_flip() + 
  geom_vline(xintercept = 1.5, linetype = "dashed",color = "red",size = .5) +
  theme_minimal() 
did_plot

# Putting together the two graphs 
did_both <- grid.arrange(
  partrends,
  did_plot,
  nrow = 1,
  ncol = 2
)
did_both

ggsave("did_both.png", plot = did_both, width = 28, height = 15, units = "cm")

# Plotting differences in sociodemographic characteristics of UKIP voters in 2010 and 2015

did_sociodem <- read.dta13("ukip_graph_sociodem.dta", encoding = 'UTF-8')

# Adding line breaks
levels(did_sociodem$Model) <- gsub(" ", "\n", levels(did_sociodem$Model))

# Ordering the variables
did_sociodem$Model <- factor(did_sociodem$Model, 
                             levels = did_sociodem$Model[order(did_sociodem$order, decreasing = TRUE)])

# Making the plot
did_sociodem_plot <- ggplot(did_sociodem, aes(x = Model, y = coef)) +   
  geom_hline(aes(yintercept = 0), linetype = 2) + 
  coord_flip(ylim = c(-1, 1)) + 
  theme_minimal(base_size = 9) +
  theme(panel.background = element_rect(fill = NA, color = "black")) +
  geom_linerange(aes(ymin = ci_lower_95, ymax = ci_upper_95), position = position_dodge(width = 0.75), size = 0.5) + 
  geom_linerange(aes(ymin = ci_lower_90, ymax = ci_upper_90), position = position_dodge(width = 0.75), size = 1) + 
  geom_point(size = 1, position = position_dodge(width = 0.75) ) + 
  ylab(" ") + 
  xlab("")
did_sociodem_plot

ggsave("ukip_did_sociodem.png", plot = did_sociodem_plot, width = 18, height = 15, units = "cm")

####
## ANALYSES FOR EP ELECTIONS ##
####

# Parallel trends graph
# Importing dataset
parallel_ep <- read.dta13("partrends_ep.dta", encoding = 'UTF-8')

# Making sure all numeric variables are such
parallel_ep$evertreat <- as.factor(parallel_ep$evertreat)
parallel_ep$time_period <- as.numeric(parallel_ep$time_period)

Palette2 <- c("blue","black")

partrends_ep = ggplot(parallel_ep, aes(x = time_period, y = meanrepvote, colour = evertreat)) +  
  geom_line(aes(linetype = evertreat, color = evertreat), position = position_dodge(0.25)) + 
  geom_point(position = position_dodge(0.25)) +
  scale_y_continuous(limits = c(0, 1.7), breaks = c(0, 0.25, 0.5, 0.75, 1, 1.25, 1.5))  + 
  ylab("Proportion of the official vote reported in the EES") +
  scale_color_manual(name = "Type of Party",values = Palette2, labels = c( ), guide = "none") +
  scale_linetype_manual(name = " ", values = c(1,2), guide = "none")+
  scale_x_continuous(name = "Election", breaks = c(1, 2, 3, 4), 
                     labels = c("1999", "2004", "2009", "2014")) +
  annotate("text", x = 2, y = 0.55, label = "Control: Conservative, Labour and LibDem", size = 3) +
  annotate("text", x = 2.3, y = 1.2, label = "Treated: UKIP", size = 3) +
  geom_errorbar(aes(ymin = ci_lower, ymax = ci_upper), width = 0.1, position = position_dodge(0.25)) +  
  labs(fill = " ") +
  theme_classic() + 
  geom_vline(xintercept = 3.85, linetype = "dashed",color = "grey",size = .5)
ggtitle("Trends in the control and treatment group") 
partrends_ep

# Plotting the findings of the DiD model

did_ep <- read.dta13("ukip_graph_ep.dta", encoding = 'UTF-8')

# Adding line breaks
levels(did_ep$Model) <- gsub("With demographic", "With\ndemographic", levels(did_ep$Model))
levels(did_ep$Model) <- gsub("demographic weights", "demographic\nweights", levels(did_ep$Model))

# Ordering the variables
did_ep$Model <- factor(did_ep$Model, 
                       levels = did_ep$Model[order(did_ep$order, decreasing = TRUE)])

did_ep_plot <- did_ep %>%
  mutate(Model = fct_reorder(Model, desc(order))) %>%
  ggplot (aes(x = Model, y = coef)) +
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +
  geom_linerange(aes(ymin = ci_lower_95, ymax = ci_upper_95), position = position_dodge(width = 0.75), size = 0.5) + 
  geom_linerange(aes(ymin = ci_lower_90, ymax = ci_upper_90), position = position_dodge(width = 0.75), size = 1) +
  geom_point(size = 1, position = position_dodge(width = 0.75) ) + 
  scale_y_continuous(name = " ", limits = c(-0.6, 0.6), breaks = c(-0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6)) +
  coord_flip() + 
  xlab(" ")  +
  geom_vline(xintercept = 1.5, linetype = "dashed",color = "red",size = .5) +
  theme_minimal() +
  ggtitle("Difference-in-differences estimates")
did_ep_plot

# Putting together the two graphs 
did_both_ep <- grid.arrange(
  partrends_ep,
  did_ep_plot,
  nrow = 1,
  ncol = 2
)
did_both_ep

ggsave("did_both_ep.png", plot = did_both_ep, width = 28, height = 15, units = "cm")

## Plotting results with different control groups

controls <- read.dta13("ukip_controls.dta")

did_controls_plot <- ggplot(controls) +
  geom_histogram(aes (x = coef), bins = 15, fill = "steelblue") +
  scale_x_continuous(name = "Effect size", limits = c(0, 0.5)) +
  scale_y_continuous(name = "Count") +
  theme_bw() +
  geom_density(aes (x = coef), fill = "#FFDB6D", color = "#C4961A", alpha = 0.5) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "red")
did_controls_plot

ggsave("did_controls.png", plot = did_controls_plot, width = 20, height = 15, units = "cm")

# Looking at the median and mean effect size
mean(controls$coef)
median(controls$coef)

## PARALLEL TRENDS PLOTS USING DIFFERENT CONTROL GROUPS
# LibDem
rawdata <- read.dta13("study3_generalelec.dta")

libdem_control <- rawdata %>%
  filter(party_code %in% c(1, 4)) %>%
  mutate(Party = ifelse(evertreat == 1, "UKIP", "LibDem")) %>%
  ggplot(aes(x = year, y = reported_vote, color = Party)) +
  geom_point() +
  scale_color_manual(values = c("blue", "red")) +
  stat_summary(fun = mean, geom = "point", position = position_dodge(width = 0)) +
  stat_summary(fun = mean, geom = "line", position = position_dodge(width = 0)) +
  scale_y_continuous(limits = c(0, 1.3)) +
  scale_x_continuous(breaks = c(2005, 2010, 2015)) +
  theme_classic() + 
  geom_vline(xintercept = 2014.5, linetype = "dashed",color = "grey",size = .5) +
  ggtitle("LibDem as control group") +
  labs(x = 'Election', y = 'Proportion of the official vote reported in post-electoral surveys') 
libdem_control

# Conservatives

conservatives_control <- rawdata %>%
  filter(party_code %in% c(1, 3)) %>%
  mutate(Party = ifelse(evertreat == 1, "UKIP", "Conservatives")) %>%
  ggplot(aes(x = year, y = reported_vote, color = Party)) +
  geom_point() +
  scale_color_manual(values = c("black", "red")) +
  stat_summary(fun = mean, geom = "point", position = position_dodge(width = 0)) +
  stat_summary(fun = mean, geom = "line", position = position_dodge(width = 0)) +
  scale_y_continuous(limits = c(0, 1.3)) +
  scale_x_continuous(breaks = c(2005, 2010, 2015)) +
  theme_classic() + 
  geom_vline(xintercept = 2014.5, linetype = "dashed",color = "grey",size = .5) +
  ggtitle("Conservatives as control group") +
  labs(x = 'Election', y = 'Proportion of the official vote reported in post-electoral surveys') 
conservatives_control

# Labour

labour_control <- rawdata %>%
  filter(party_code %in% c(1, 2)) %>%
  mutate(Party = ifelse(evertreat == 1, "UKIP", "Labour")) %>%
  ggplot(aes(x = year, y = reported_vote, color = Party)) +
  geom_point() +
  scale_color_manual(values = c("darkgreen", "red")) +
  stat_summary(fun = mean, geom = "point", position = position_dodge(width = 0)) +
  stat_summary(fun = mean, geom = "line", position = position_dodge(width = 0)) +
  scale_y_continuous(limits = c(0, 1.3)) +
  scale_x_continuous(breaks = c(2005, 2010, 2015)) +
  theme_classic() + 
  geom_vline(xintercept = 2014.5, linetype = "dashed",color = "grey",size = .5) +
  ggtitle("Labour as control group") +
  labs(x = 'Election', y = 'Proportion of the official vote reported in post-electoral surveys') 
labour_control

# Putting the three plots together
controls <- grid.arrange(
  libdem_control,
  conservatives_control,
  labour_control,
  nrow = 1,
  ncol = 3)
controls

ggsave("ukip_controls.png", plot = controls, width = 30, height = 15, units = "cm")

